About the project

Tällä kurssilla käytetään ohjelmistoja Rstudio, Github ja Github desktop. Rstudiota käytetään tilastolliseen analyysiin ja Githubia tiedostojen säilömiseen. Github desktopilla yhdistetään tietokoneen tiedostot ja säilytyspaikka Githubissa.

https://github.com/aolkinuo/IODS-project


Regression and model validation

learning2014=read.table("data/learning2014.txt", sep = ",", header = TRUE)

str(learning2014)
## 'data.frame':    166 obs. of  7 variables:
##  $ gender  : Factor w/ 2 levels "F","M": 1 2 1 2 2 1 2 1 2 1 ...
##  $ Age     : int  53 55 49 53 49 38 50 37 37 42 ...
##  $ Attitude: int  37 31 25 35 37 38 35 29 38 21 ...
##  $ deep    : num  3.58 2.92 3.5 3.5 3.67 ...
##  $ stra    : num  3.38 2.75 3.62 3.12 3.62 ...
##  $ surf    : num  2.58 3.17 2.25 2.25 2.83 ...
##  $ Points  : int  25 12 24 10 22 21 21 31 24 26 ...
dim(learning2014)
## [1] 166   7

Aineistossa on 166 havaintoyksikköä ja seitsemän muuttujaa. Muuttujat Age, Attitude, Points ovat numeerisia kokonaislukumuuttujia ja ne kertovat henkilön iän, asenteen tilastotiedettä kohtaan ja koepisteet. Ikä on laskettu syntymävuoden avulla. Muuttujat deep, stra ja surf ovat numeerisia muuttujia. Ne kertovat henkilön vastaukset syvällisiin, strategisiin ja pinnallisiin kysymyksiin. gender on luokiteltu muuttuja ja sen arvo 1 kuvaa naista ja arvo 2 miestä.

summary(learning2014)
##  gender       Age           Attitude          deep            stra      
##  F:110   Min.   :17.00   Min.   :14.00   Min.   :1.583   Min.   :1.250  
##  M: 56   1st Qu.:21.00   1st Qu.:26.00   1st Qu.:3.333   1st Qu.:2.625  
##          Median :22.00   Median :32.00   Median :3.667   Median :3.188  
##          Mean   :25.51   Mean   :31.43   Mean   :3.680   Mean   :3.121  
##          3rd Qu.:27.00   3rd Qu.:37.00   3rd Qu.:4.083   3rd Qu.:3.625  
##          Max.   :55.00   Max.   :50.00   Max.   :4.917   Max.   :5.000  
##       surf           Points     
##  Min.   :1.583   Min.   : 7.00  
##  1st Qu.:2.417   1st Qu.:19.00  
##  Median :2.833   Median :23.00  
##  Mean   :2.787   Mean   :22.72  
##  3rd Qu.:3.167   3rd Qu.:27.75  
##  Max.   :4.333   Max.   :33.00

Aineiston yhteenvedosta, joka on tämän tekstin yläpuolella nähdään muuttujien Age, Attitude, deep, stra, surf ja Points minimit, maksimit, ensimmäinen ja kolmas kvantiili, mediaani ja keskiarvo. Muuttujasta gender nähdään sen kahden arvon frekvenssit.

library(ggplot2)

library(GGally)

pairs(learning2014[-1], col = learning2014$gender)

p=ggpairs(learning2014, mapping = aes(col = gender, alpha = 0.3), lower = list(combo = wrap("facethist", bins = 20)))

p

Kuvasta nähdään että muuttujat Attitude, deep ja stra ovat jakautuneet lähes normaalisti naisten ja miesten joukossa. Myös surf on naisten osalta jakautunut lähes normaalisti. Age on vinoutunut, skewed, oikealle ja myös surf on vinoutunut oikealle miesten osalta. Points-muuttujan jakauma on vinoutunut vasemmalle. Korrelaatiokertoimista nähdään että muuttujat Attitude ja Age, deep ja Age, deep ja Attitude, stra ja Attitude, Points ja Attitude, stra ja deep, Points ja deep korreloivat suhteellisen paljon naisten joukossa. Muuttujat Attitude ja Age, deep ja Age, stra ja Age, surf ja Attitude, Points ja Attitude, surf ja deep, Points ja deep korreloivat suhteellisen paljon miesten joukossa. Suuren korrelaatiokertoimen tulkinta on esimerkiksi muuttujien Attitude ja Age kohdalla se että asenteen tilastotiedettä kohtaan ja iän välillä on riippuvuutta. Jos korrelaatiokerroin olisi yksi, toisen muuttujan arvon voisi laskea tarkasti, lineaarisesti toisen arvosta.

malli1=lm(Points~Age+Attitude+deep, data=learning2014)

summary(malli1)
## 
## Call:
## lm(formula = Points ~ Age + Attitude + deep, data = learning2014)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.206  -3.430   0.204   3.979  10.950 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 15.60773    3.38966   4.605 8.32e-06 ***
## Age         -0.07716    0.05322  -1.450    0.149    
## Attitude     0.35941    0.05696   6.310 2.56e-09 ***
## deep        -0.60275    0.75031  -0.803    0.423    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.307 on 162 degrees of freedom
## Multiple R-squared:  0.2043, Adjusted R-squared:  0.1896 
## F-statistic: 13.87 on 3 and 162 DF,  p-value: 4.305e-08

Kun estimoidaan pienimmän neliösumman menetelmällä malli jossa ikä, asenne tilastotiedettä kohtaan ja syvälliset kysymykset selittävät koepisteitä, saadaan tulokseksi tulostus joka on tämän tekstin yläpuolella.

Ensimmäisessä sarakkeessa jossa on numeroita näkyy selittäviin muuttujiin liittyvien parametrien eli parametrivektorin arvot. Niistä nähdään että kun ikä kasvaa yhdellä, koepisteet laskevat noin -0.08 verran. Kun asenne tilastotiedettä kohtaan nousee yhdellä, koepisteet kasvavat noin 0.36 verran. Kun taas syvälliset kysymykset nousevat yhdellä, koepisteet laskevat 0.6 verran. Kun ikä, asenne ja syvälliset kysymykset ovat nolla, koepisteet saavan vakiotermin mukaisen arvon eli noin 15.61.

Tähdet kuvassa kertovat selittävien muuttujien tilastollisesta merkitsevyydestä. Tähtien taustalla on tilastollinen testi jossa testataan sitä voidaanko nollahypoteesi siitä että selittävään muuttujaan liittyvä parametri on nolla hylätä. Nollahypoteesit että vakiotermin ja asennemuuttujan parametrit ovat nolla voidaan hylätä 0.1 prosentin merkitsevyystasolla eli asennemuuttujan ja vakiotermin parametrit eli vaikutukset selitettävään koepistemuuttujaan ovat merkitseviä.Nollahypoteeseja siitä että asenne- ja syvälliset kysymykset-muuttujien parametrit ovat nolla, ei voida hylätä. Asenne- ja syvälliset kysymykset -muuttujien vaikutukset koepistemuuttujaan eivät siis ole merkitseviä. Poistetaan muuttujien vaikutusten huonon merkitsevyyden puolesta asenne- ja syvälliset #kysymykset-muuttujat mallista ja estimoidaan pienimmän neliösumman #menetelmällä malli jossa on yksi selittävä muuttuja. Se on asenne tilastotiedettä kohtaan.

malli2=lm(Points~Attitude, data=learning2014)

summary(malli2)
## 
## Call:
## lm(formula = Points ~ Attitude, data = learning2014)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.9763  -3.2119   0.4339   4.1534  10.6645 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.63715    1.83035   6.358 1.95e-09 ***
## Attitude     0.35255    0.05674   6.214 4.12e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.32 on 164 degrees of freedom
## Multiple R-squared:  0.1906, Adjusted R-squared:  0.1856 
## F-statistic: 38.61 on 1 and 164 DF,  p-value: 4.119e-09

Tulostuksesta nähdään että asennemuuttujan parametrin arvo on noin 0.35. Se tarkoittaa että kun asennemuuttuja kasvaa yhdellä, koepisteet kasvavat noin 0.35 verran. Asennemuuttujaan ja vakiotermiin liittyvät parametrit ovat tilastollisesti merkitseviä 0.1 prosentin merkitsevyystasoilla.

Multiple R-squared ja Adjusted R-squared ovat mallin selitysasteita. Selitysaste kertoo sen kuinka paljon malli, eli selittävät muuttujat jotka on mallinnettu käyttämällä lineaarista regressiota, selittää selitettävän muuttujan vaihtelusta. Yläpuolella olevan mallin selitysaste kertoo siis sen että asenteen tilastotiedettä kohtaan ja vakiotermin vaihtelut selittävät noin 19 prosenttia koepisteiden vaihtelusta. Selitysaste kertoo tavallaan mallin suorituskyvystä ja se saa arvon nollan ja yhden välillä.

Adjusted R-squared saa aina pienemmän arvon kuin Multiple R-squared sillä siinä huomioidaan mallin monimutkaisuus eli muuttujien määrä. Adjusted R-squared on tarkempi mallin suorituskyvyn mittari kuin Multiple R-squared. Uuden selittävän muuttujan lisääminen malliin kasvattaa todennäköisesti Multiple R-squaredia ja pienentää Adjusted R-squaredia.

par(mfrow = c(2,2))
plot(malli2, which = c(1,2,5))

On piirretty kolme kuvaa, kuva jossa on mallin kaksi sovite ja residuaalit, QQ-plot-kuva jossa on teoreettiset kvantiilit ja standardoidut residuaalit ja kuva jossa on vipuvoima ja standardoidut residuaalit.

Ensimmäisestä kuvasta eli kuvasta jossa on sovite x-akselilla ja residuaalit y-akselilla voidaan tutkia sitä onko virheillä vakio varianssi. Koska havainnot ovat jakautuneet suhteellisen tasaisesti nollasuoran ympärille eivätkä muodosta jotakin muotoa tai mallia, voidaan sanoa että virhetermin varianssi on vakio eikä virheiden suuruus riipu selittävistä muuttujista. Virheiden vakio varianssi on yksi oletus joka usein tehdään kun estimoidaan lineaarinen regressiomalli. Se näyttäisi mallin kaksi tapauksessa olevan voimassa.

QQ-plot-kuvasta eli kuvasta jossa on teoreettiset kvantiilit x-akselilla ja standardoidut residuaalit y-akselilla voidaan tehdä päätelmiä virheiden normaalisuudesta. Koska havainnot ovat QQ-plot-kuvassa suoran päällä, virhetermin voidaan sanoa noudattavan suurin piirtein normaalijakaumaa. Niin usein halutaan olettaa ja nyt mallin kaksi tapauksessa voidaan sanoa että oletus että virheet ovat normaalisti jakaantuneet voidaan pitää voimassa.

Kolmannesta kuvasta, jossa on vipuvoima eli Cookin etäisyys x-akselilla ja standardoidut residuaalit y-akselilla voidaan päätellä vipuvoiman suuruus. Vipuvoima mittaa sitä kuinka paljon yhdellä havainnolla on vaikutusta malliin. Hajontakuva jossa on residuaalit ja vipuvoima voi auttaa tunnistamaan ne havainnot joilla on poikkeuksellisen paljon vaikutusta malliin. Kuvasta nähdään että havainto joka saa suurimman vipuvoiman tai Cookin etäisyyden arvon saa suunnilleen arvon 0.05. Luku on sen verran pieni että vipuvoiman mallissa kaksi voidaan sanoa olevan tavallista tai pientä.

Muita lineaarisen regressiomallin oletuksia ovat virheiden normaalisuuden ja vakion varianssin lisäksi oletukset että virheet eivät korreloi eivätkä riipu selittävistä muuttujista.


alc=read.table("data/alc.txt", sep = ",", header = TRUE)

head(alc)
##   school sex age address famsize Pstatus Medu Fedu     Mjob     Fjob
## 1     GP   F  18       U     GT3       A    4    4  at_home  teacher
## 2     GP   F  17       U     GT3       T    1    1  at_home    other
## 3     GP   F  15       U     LE3       T    1    1  at_home    other
## 4     GP   F  15       U     GT3       T    4    2   health services
## 5     GP   F  16       U     GT3       T    3    3    other    other
## 6     GP   M  16       U     LE3       T    4    3 services    other
##       reason nursery internet guardian traveltime studytime failures
## 1     course     yes       no   mother          2         2        0
## 2     course      no      yes   father          1         2        0
## 3      other     yes      yes   mother          1         2        2
## 4       home     yes      yes   mother          1         3        0
## 5       home     yes       no   father          1         2        0
## 6 reputation     yes      yes   mother          1         2        0
##   schoolsup famsup paid activities higher romantic famrel freetime goout
## 1       yes     no   no         no    yes       no      4        3     4
## 2        no    yes   no         no    yes       no      5        3     3
## 3       yes     no  yes         no    yes       no      4        3     2
## 4        no    yes  yes        yes    yes      yes      3        2     2
## 5        no    yes  yes         no    yes       no      4        3     2
## 6        no    yes  yes        yes    yes       no      5        4     2
##   Dalc Walc health absences G1 G2 G3 alc_use high_use
## 1    1    1      3        5  2  8  8     1.0    FALSE
## 2    1    1      3        3  7  8  8     1.0    FALSE
## 3    2    3      3        8 10 10 11     2.5     TRUE
## 4    1    1      5        1 14 14 14     1.0    FALSE
## 5    1    2      5        2  8 12 12     1.5    FALSE
## 6    1    2      5        8 14 14 14     1.5    FALSE
attach(alc)

Aineiston muuttujia ovat esimerkiksi alkoholin käyttö viikolla, alkoholin käyttö viikonloppuna, terveys, poissaolot, alkoholin käyttö koko viikon aikana, suuri alkoholinkulutus, opiskeluun käytetty aika, huoltaja, perheen koko, osoite, ikä, sukupuoli ja koulu. Siinä oli vain osa muuttujista. Loput näkyvät tulosteesta.

Tutkin seuraavaksi suuren alkoholinkäytön joka on looginen muuttuja yhteyttä terveyteen, poissaoloihin, opiskeluun käytettyyn aikaan ja sukupuoleen. Hypoteesini on että terveys ja suuri alkoholinkäyttö korreloivat negatiivisesti. Arvioin myös että opiskeluun käytetyllä ajalla ja suurella alkoholinkäytöllä on negatiivinen korrelaatio. Arvioin myös että poissaolot ja suuri alkoholinkäyttö korreloivat positiivisesti. Hypoteesini on että sukupuoli vaikuttaa suureen alkoholinkäyttöön niin että miehissä on enemmän heitä jotka käyttävät paljon alkoholia eli heitä joille high_use-muuttuja saa arvon TRUE.

table(health,high_use)
##       high_use
## health FALSE TRUE
##      1    35   11
##      2    28   16
##      3    62   19
##      4    49   18
##      5    94   50
table(absences,high_use)
##         high_use
## absences FALSE TRUE
##       0     52   13
##       1     38   13
##       2     42   16
##       3     33    8
##       4     24   12
##       5     16    6
##       6     16    5
##       7      9    3
##       8     14    6
##       9      6    6
##       10     5    2
##       11     2    4
##       12     4    4
##       13     1    1
##       14     1    6
##       16     0    1
##       17     0    1
##       18     1    1
##       19     0    1
##       20     2    0
##       21     1    1
##       26     0    1
##       27     0    1
##       29     0    1
##       44     0    1
##       45     1    0
table(studytime,high_use)
##          high_use
## studytime FALSE TRUE
##         1    58   42
##         2   135   60
##         3    52    8
##         4    23    4
table(sex,high_use)
##    high_use
## sex FALSE TRUE
##   F   156   42
##   M   112   72
barplot(health,high_use)

barplot(absences,high_use)

barplot(studytime,high_use)

boxplot(health,high_use)

boxplot(absences,high_use)

boxplot(studytime,high_use)

boxplot(sex,high_use)

Laatikkokuvista nähdään että hypoteesini terveyden ja suuren alkoholinkäytön ja opiskeluun käytetyn ajan ja suuren alkoholinkäytön välisistä suhteista olivat oikein. Niiden joiden alkoholinkäyttö ei ole suurta terveyden mediaani ja ala- ja yläkvantiilit ovat selvästi korkeammalla tasolla kuin heidän joiden alkoholinkäyttö on suurta vastaavat tunnusluvut. Samoin niiden joiden alkoholinkäyttö ei ole suurta opiskeluun käytetyn ajan mediaani ja ala- ja yläkvantiilit ovat korkeammalla tasolla kuin heidän joiden alkoholinkäyttö on suurta vastaavat tunnusluvut. Alkoholinkäytön kasvaminen siis vaikuttaa terveyteen ja opiskeluun käytettyyn aikaan negatiivisesti. Poissaolojen suhteen arvioni meni pieleen. Heidän joiden alkoholinkäyttö ei ole suurta poissaolojen mediaani ja ylä- ja alakvantiilit ovat korkeammalla tasolla kuin heidän joiden alkoholinkäyttö on suurta vastaavat tunnusluvut. Alkoholinkäytön kasvaminen siis vaikuttaa poissaoloihin negatiivisesti mikä oli päinvastoin kuin mitä arvioin. Ristiintaulukosta nähdään että sellaisia miehiä joiden alkoholinkäyttö on suurta on enemmän kuin sellaisia naisia joiden alkoholinkäyttö on suurta. Sellaisia miehiä joiden alkoholinkäyttö ei ole suurta on vähemmän kuin sellaisia naisia joiden alkoholinkäyttö ei ole suurta. Miehissä on siis oltava suhteellisesti suurempi osuus heitä joiden alkoholinkäyttö on suurta kuin naisissa. Hypoteesini oli sukupuolen ja suuren alkoholinkäytön suhteen osalta oikea.

m = glm(high_use ~ health + absences + studytime + sex, data = alc, family = "binomial")

summary(m)
## 
## Call:
## glm(formula = high_use ~ health + absences + studytime + sex, 
##     family = "binomial", data = alc)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2074  -0.8367  -0.5983   1.0910   2.1407  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.08010    0.52055  -2.075 0.037995 *  
## health       0.04493    0.08631   0.521 0.602675    
## absences     0.08881    0.02307   3.850 0.000118 ***
## studytime   -0.39817    0.15960  -2.495 0.012603 *  
## sexM         0.78234    0.25261   3.097 0.001954 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 465.68  on 381  degrees of freedom
## Residual deviance: 423.00  on 377  degrees of freedom
## AIC: 433
## 
## Number of Fisher Scoring iterations: 4
coef(m)
## (Intercept)      health    absences   studytime        sexM 
## -1.08009806  0.04492699  0.08881324 -0.39817186  0.78234486

Logistisen regressiomallin estimoinnin tulostuksesta nähdään että kun terveys kasvaa yhdellä, logaritmi suuren alkoholinkäytön todennäköisyyden ja yksi miinus suuren alkoholinkäytön todennäköisyyden osamäärästä eli log(p/(1-p)):stä kasvaa 0,04 verran. Samoin jos poissaolot kasvavat yhdellä, log(p/(1-p)) kasvaa noin 0,09 verran. Opiskeluun käytetty aika -muuttujan kerroin tulkitaan vastaavalla tavalla. Sukupuoli on luokiteltu muuttuja ja sen kertoimen tulkinta on erilainen kuin muiden muuttujien kertoimien tulkinta. R on valinnut naiset vertailuryhmäksi ja miehet ovat selittävä muuttuja. Vertailuryhmän eli naisten kerroin on mallin vakiotermi. Miesten kerroin kuvaa eroa vakiotermiin eli naisten kertoimeen. Miesten todellinen kerroin on vakiotermi-miesten kerroin eli noin -1.08-0.78=-0.3. Selittäjistä tilastollisesti merkitsevä kertoimen estimaatti on poissaoloilla, opiskeluun käytetyllä ajalla, miehillä ja vakiotermillä. Se että miesten kertoimen estimaatti on tilastollisesti merkitsevä tarkoittaa että hypoteesi että vakiotermin eli vertailuryhmän kertoimen ja miesten kertoimen erotus on nolla voidaan hylätä yhden prosentin merkitsevyystasolla. Terveyden kertoimen estimaatti ei ole tilastollisesti merkitsevä eli hypoteesia siitä että terveyden kerroin on nolla ei voida hylätä yleisillä merkitsevyystasoilla.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:GGally':
## 
##     nasa
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
OR =coef(m) %>% exp

CI=exp(confint(m))
## Waiting for profiling to be done...
cbind(OR, CI)
##                    OR     2.5 %    97.5 %
## (Intercept) 0.3395622 0.1203826 0.9317189
## health      1.0459515 0.8842249 1.2412355
## absences    1.0928765 1.0468090 1.1460148
## studytime   0.6715466 0.4866129 0.9113044
## sexM        2.1865935 1.3377030 3.6085532

Terveyden odds ratio kertoo suuren alkoholinkäytön todennäköisyyden henkilöillä joilla terveys-muuttuja muuttuu yhdellä yksiköllä ja suuren alkoholinkäytön todennäköisyyden henkilöillä jolla terveys-muuttuja ei muutu yhdellä yksiköllä osamäärän. Koska se on suurempaa kuin yksi, terveyden muutos yhdellä yksiköllä liittyy positiivisesti suureen alkoholinkäyttöön. Tämä on eri tulos kuin hypoteesissani arvioin. Myös poissaolojen odds ratio on yli yksi eli poissaolojen muutos yhdellä yksiköllä liittyy positiivisesti suureen alkoholinkäyttöön. Se on hypoteesini mukainen tulos. Opiskeluun käytetyn ajan odds ratio on alle yksi mikä tarkoittaa että opiskeluun käytetyn ajan muutos yhdellä yksiköllä liittyy negatiivisesti suureen alkoholinkäyttöön. Se on arvioni mukainen tulos. Se että on mies näyttää liittyvän positiivisesti suureen alkoholinkäyttöön, sillä miesten odds ratio on yli kaksi. Se on arvioni mukaista.

#Lasketaan malli ilman terveysmuuttujaa jonka kertoimen estimaatti ei ollut tilastollisesti merkitsevä. 
m = glm(high_use ~ absences + studytime + sex, data = alc, family = 
"binomial")

#Ennustetaan suuren alkoholinkäytön todennäköisyys.
probabilities = predict(m, type = "response")

#Lisätään ennustetut todennäköisyydet aineistoon.
alc = mutate(alc, probability = probabilities)

#Käytetään todennäköisyyksiä suuren alkoholinkäytön ennustamiseen,
alc = mutate(alc, prediction = probabilities>0.5)

#Ristiintaulukoidaan korkea alkoholinkäyttö-muuttuja ennusteiden muuttujasta kanssa.
table(high_use = alc$high_use, prediction = probabilities>0.5)
##         prediction
## high_use FALSE TRUE
##    FALSE   258   10
##    TRUE     88   26

Ristiintaulukosta nähdään että suuri alkoholinkäyttö-muuttujan arvot poikkeavat siitä logit-mallin avulla tehdyistä ennusteista vähän. 10 kertaa ennuste on suuri alkoholinkäyttö silloin kun muuttuja saa ei suurta alkoholinkäyttöä vastaavan arvon. Samoin 88 kertaa ennuste on ei suurta alkoholinkäyttöä kun muuttuja saa suurta alkoholinkäyttöä vastaavan arvon.

# Haetaan paketit dplyr ja ggplot2.
library(dplyr); library(ggplot2)

# Piirretään hajontakuva korkeasta alkoholin käytöstä ja sen ennustetuista todennäköisyyksistä.
g = ggplot(alc, aes(x = probability, y = high_use, col=prediction))

# Määritellään geom pisteinä ja piirretään uusi kuva. 
g+geom_point()

#Taulukoidaan suuri alkoholinkäyttö ja ennusteet siitä. 
table(high_use = alc$high_use, prediction = alc$prediction)%>%prop.table()%>%addmargins()
##         prediction
## high_use      FALSE       TRUE        Sum
##    FALSE 0.67539267 0.02617801 0.70157068
##    TRUE  0.23036649 0.06806283 0.29842932
##    Sum   0.90575916 0.09424084 1.00000000
# Määritellään loss function eli keskimääräinen ennustevirhe. 
loss_func = function(class, prob) {
  n_wrong = abs(class - prob) > 0.5
  mean(n_wrong)
}

# Kutsutaan loss funktio jotta voidaan laskea väärien ennusteiden keskimääräinen määrä aineistossa. 
loss_func(class = alc$high_use, prob = alc$probability)
## [1] 0.2565445

Väärin luokiteltuja henkilöitä on noin 26 prosenttia kaikista henkilöistä. Ylemmästä taulukosta nähdään että ennusteen mukaan henkilöitä joilla ei ole suurta alkoholinkäyttöä on noin 91 prosenttia kaikista henkilöistä. Muuttujan arvojen mukaan henkilöitä joilla ei ole suurta alkoholinkäyttöä on noin 70 prosenttia kaikista henkilöistä. Huomataan taas että ennuste eroaa muuttujan arvoista.

``` ***

# access the MASS package
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
# load the data
data("Boston")

# explore the dataset
str(Boston)
## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
dim(Boston)
## [1] 506  14

We loaded the Boston data from the MASS package. We explored the structure and the dimensions of the data. From the output of the r command dim we see that the data contains 14 variables and 506 observations. From the output of r command str we see that all variables are numerical and some of them get integer values. From the web page concerning the Boston dataset we saw that the variables are the following

crim

per capita crime rate by town.

zn

proportion of residential land zoned for lots over 25,000 sq.ft.

indus

proportion of non-retail business acres per town.

chas

Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).

nox

nitrogen oxides concentration (parts per 10 million).

rm

average number of rooms per dwelling.

age

proportion of owner-occupied units built prior to 1940.

dis

weighted mean of distances to five Boston employment centres.

rad

index of accessibility to radial highways.

tax

full-value property-tax rate per \$10,000.

ptratio

pupil-teacher ratio by town.

black

1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town.

lstat

lower status of the population (percent).

medv

median value of owner-occupied homes in \$1000s.
summary(Boston)
##       crim                zn             indus            chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00
library(ggplot2)

library(GGally)

pairs(Boston[-1])

p=ggpairs(Boston, mapping = aes(alpha = 0.3), lower = list(combo = wrap("facethist", bins = 20)))

p

# calculate the correlation matrix
cor_matrix<-cor(Boston) 

# print the correlation matrix
cor_matrix
##                crim          zn       indus         chas         nox
## crim     1.00000000 -0.20046922  0.40658341 -0.055891582  0.42097171
## zn      -0.20046922  1.00000000 -0.53382819 -0.042696719 -0.51660371
## indus    0.40658341 -0.53382819  1.00000000  0.062938027  0.76365145
## chas    -0.05589158 -0.04269672  0.06293803  1.000000000  0.09120281
## nox      0.42097171 -0.51660371  0.76365145  0.091202807  1.00000000
## rm      -0.21924670  0.31199059 -0.39167585  0.091251225 -0.30218819
## age      0.35273425 -0.56953734  0.64477851  0.086517774  0.73147010
## dis     -0.37967009  0.66440822 -0.70802699 -0.099175780 -0.76923011
## rad      0.62550515 -0.31194783  0.59512927 -0.007368241  0.61144056
## tax      0.58276431 -0.31456332  0.72076018 -0.035586518  0.66802320
## ptratio  0.28994558 -0.39167855  0.38324756 -0.121515174  0.18893268
## black   -0.38506394  0.17552032 -0.35697654  0.048788485 -0.38005064
## lstat    0.45562148 -0.41299457  0.60379972 -0.053929298  0.59087892
## medv    -0.38830461  0.36044534 -0.48372516  0.175260177 -0.42732077
##                  rm         age         dis          rad         tax
## crim    -0.21924670  0.35273425 -0.37967009  0.625505145  0.58276431
## zn       0.31199059 -0.56953734  0.66440822 -0.311947826 -0.31456332
## indus   -0.39167585  0.64477851 -0.70802699  0.595129275  0.72076018
## chas     0.09125123  0.08651777 -0.09917578 -0.007368241 -0.03558652
## nox     -0.30218819  0.73147010 -0.76923011  0.611440563  0.66802320
## rm       1.00000000 -0.24026493  0.20524621 -0.209846668 -0.29204783
## age     -0.24026493  1.00000000 -0.74788054  0.456022452  0.50645559
## dis      0.20524621 -0.74788054  1.00000000 -0.494587930 -0.53443158
## rad     -0.20984667  0.45602245 -0.49458793  1.000000000  0.91022819
## tax     -0.29204783  0.50645559 -0.53443158  0.910228189  1.00000000
## ptratio -0.35550149  0.26151501 -0.23247054  0.464741179  0.46085304
## black    0.12806864 -0.27353398  0.29151167 -0.444412816 -0.44180801
## lstat   -0.61380827  0.60233853 -0.49699583  0.488676335  0.54399341
## medv     0.69535995 -0.37695457  0.24992873 -0.381626231 -0.46853593
##            ptratio       black      lstat       medv
## crim     0.2899456 -0.38506394  0.4556215 -0.3883046
## zn      -0.3916785  0.17552032 -0.4129946  0.3604453
## indus    0.3832476 -0.35697654  0.6037997 -0.4837252
## chas    -0.1215152  0.04878848 -0.0539293  0.1752602
## nox      0.1889327 -0.38005064  0.5908789 -0.4273208
## rm      -0.3555015  0.12806864 -0.6138083  0.6953599
## age      0.2615150 -0.27353398  0.6023385 -0.3769546
## dis     -0.2324705  0.29151167 -0.4969958  0.2499287
## rad      0.4647412 -0.44441282  0.4886763 -0.3816262
## tax      0.4608530 -0.44180801  0.5439934 -0.4685359
## ptratio  1.0000000 -0.17738330  0.3740443 -0.5077867
## black   -0.1773833  1.00000000 -0.3660869  0.3334608
## lstat    0.3740443 -0.36608690  1.0000000 -0.7376627
## medv    -0.5077867  0.33346082 -0.7376627  1.0000000

We showed a graphical overview of the data and summaries of the variables in the data. From the summaries we can see that for example minimum of crime rate in one town is about 0.006, median is about 0.257, mean is about 3.614 and maximum is about 88.976. From the graph we can see that the distributions of one variable against another are rarely normal. Only variable average number of rooms per dwelling´s distribution when the same variable, average number of rooms per dwelling, is on the y-axis looks normal. From correlation matrix we see that for example crime rate in one town correlates a lot with index of accessibility to radial highways and with full-value property-tax rate per 10,000 dollars. Correlation of crime rate with index of accessibility to radial highways is 0.63. Correlation of crime rate and full-value property-tax rate per 10,000 dollars is 0.58.

# center and standardize variables
boston_scaled <- scale(Boston)

# summaries of the scaled variables
summary(boston_scaled)
##       crim                 zn               indus        
##  Min.   :-0.419367   Min.   :-0.48724   Min.   :-1.5563  
##  1st Qu.:-0.410563   1st Qu.:-0.48724   1st Qu.:-0.8668  
##  Median :-0.390280   Median :-0.48724   Median :-0.2109  
##  Mean   : 0.000000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.007389   3rd Qu.: 0.04872   3rd Qu.: 1.0150  
##  Max.   : 9.924110   Max.   : 3.80047   Max.   : 2.4202  
##       chas              nox                rm               age         
##  Min.   :-0.2723   Min.   :-1.4644   Min.   :-3.8764   Min.   :-2.3331  
##  1st Qu.:-0.2723   1st Qu.:-0.9121   1st Qu.:-0.5681   1st Qu.:-0.8366  
##  Median :-0.2723   Median :-0.1441   Median :-0.1084   Median : 0.3171  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.:-0.2723   3rd Qu.: 0.5981   3rd Qu.: 0.4823   3rd Qu.: 0.9059  
##  Max.   : 3.6648   Max.   : 2.7296   Max.   : 3.5515   Max.   : 1.1164  
##       dis               rad               tax             ptratio       
##  Min.   :-1.2658   Min.   :-0.9819   Min.   :-1.3127   Min.   :-2.7047  
##  1st Qu.:-0.8049   1st Qu.:-0.6373   1st Qu.:-0.7668   1st Qu.:-0.4876  
##  Median :-0.2790   Median :-0.5225   Median :-0.4642   Median : 0.2746  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.6617   3rd Qu.: 1.6596   3rd Qu.: 1.5294   3rd Qu.: 0.8058  
##  Max.   : 3.9566   Max.   : 1.6596   Max.   : 1.7964   Max.   : 1.6372  
##      black             lstat              medv        
##  Min.   :-3.9033   Min.   :-1.5296   Min.   :-1.9063  
##  1st Qu.: 0.2049   1st Qu.:-0.7986   1st Qu.:-0.5989  
##  Median : 0.3808   Median :-0.1811   Median :-0.1449  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.4332   3rd Qu.: 0.6024   3rd Qu.: 0.2683  
##  Max.   : 0.4406   Max.   : 3.5453   Max.   : 2.9865
# change the object to data frame
boston_scaled=as.data.frame(boston_scaled)

We scaled the variables in the Boston dataset. From the summary of the scaled variables we see that for example statistics of crime rate have changed radically. Minimum changed from about 0.006 to about -0.419. Median changed from about 0.257 to about -0.39. Mean changed from about 3.614 to 0 and maximum changed from about 88.976 to about 9.924110.

# save the scaled crim as scaled_crim
scaled_crim <- boston_scaled$crim

# create a quantile vector of crim and print it
bins <- quantile(scaled_crim)

# create a categorical variable 'crime', label-käsky on varmaan väärin
crime <- cut(scaled_crim, breaks = bins, include.lowest = TRUE, label=c("low","med_low","med_high","high"))

# remove original crim from the dataset
boston_scaled <- dplyr::select(boston_scaled, -crim)

# add the new categorical value to scaled data
boston_scaled <- data.frame(boston_scaled, crime)

# number of rows in the Boston dataset 
n <- nrow(boston_scaled)

# choose randomly 80% of the rows
ind <- sample(n,  size = n * 0.8)

# create train set
train <- boston_scaled[ind,]

# create test set 
test <- boston_scaled[-ind,]

We have created a categorical variable of the crime rate in the Boston dataset from the scaled crime rate. We used the quantiles as the break points in the categorical variable and dropped the old crime rate variable from the dataset. We also divided the dataset to train and test sets, so that 80% of the data belongs to the train set.

# linear discriminant analysis
lda.fit <- lda(crime ~ ., data = train)

# print the lda.fit object
lda.fit
## Call:
## lda(crime ~ ., data = train)
## 
## Prior probabilities of groups:
##       low   med_low  med_high      high 
## 0.2549505 0.2425743 0.2400990 0.2623762 
## 
## Group means:
##                   zn      indus        chas        nox           rm
## low       0.98660934 -0.9375501 -0.11943197 -0.9030104  0.467705696
## med_low  -0.05496865 -0.2881305 -0.03128211 -0.5611315 -0.080692042
## med_high -0.39794966  0.2344301  0.21473488  0.4119683  0.001914263
## high     -0.48724019  1.0149946 -0.12375925  1.0540005 -0.375432529
##                 age        dis        rad        tax     ptratio
## low      -0.8950731  0.9031863 -0.6941969 -0.7372077 -0.46334026
## med_low  -0.3476199  0.4090516 -0.5412349 -0.4749902 -0.06382916
## med_high  0.4035374 -0.3969694 -0.3875099 -0.2605198 -0.26660397
## high      0.8033303 -0.8391519  1.6596029  1.5294129  0.80577843
##                black       lstat        medv
## low       0.37592101 -0.77181900  0.56847880
## med_low   0.31816897 -0.15285299  0.03138186
## med_high  0.07286901  0.02425766  0.06918111
## high     -0.72806978  0.91710547 -0.78129338
## 
## Coefficients of linear discriminants:
##                 LD1          LD2         LD3
## zn       0.11964684  0.666620434 -0.94288640
## indus   -0.01286165 -0.420516195  0.46997612
## chas    -0.06265531 -0.092893089  0.09608159
## nox      0.26323618 -0.733362331 -1.36637792
## rm      -0.07581601  0.027189675 -0.02200715
## age      0.30403875 -0.296225043 -0.15025983
## dis     -0.15537536 -0.241088177  0.37789867
## rad      3.22383311  0.822847267  0.06381472
## tax     -0.00562437  0.152887149  0.30659036
## ptratio  0.11632448  0.050901219 -0.17709758
## black   -0.14576096  0.002426148  0.10598289
## lstat    0.12536363 -0.033984115  0.55799190
## medv     0.07859804 -0.273197730 -0.11515312
## 
## Proportion of trace:
##    LD1    LD2    LD3 
## 0.9516 0.0374 0.0110
# the function for lda biplot arrows
lda.arrows <- function(x, myscale = 1, arrow_heads = 0.1, color = "orange", tex = 0.75, choices = c(1,2)){
  heads <- coef(x)
  arrows(x0 = 0, y0 = 0, 
         x1 = myscale * heads[,choices[1]], 
         y1 = myscale * heads[,choices[2]], col=color, length = arrow_heads)
  text(myscale * heads[,choices], labels = row.names(heads), 
       cex = tex, col=color, pos=3)
}

# target classes as numeric
classes <- as.numeric(train$crime)

# plot the lda results
plot(lda.fit, dimen = 2, col = classes, pch = classes)
lda.arrows(lda.fit, myscale = 1)

We fitted the linear discriminant analysis that is LDA on the train set. We used the categorical crime rate as the target variable and all the other variables in the dataset as predictor variables. We drew the LDA (bi)plot.

# save the correct classes from test data
correct_classes <- test$crime

# remove the crime variable from test data
test <- dplyr::select(test, -crime)

# predict classes with test data
lda.pred <- predict(lda.fit, newdata = test)

# cross tabulate the results
table(correct = correct_classes, predicted = lda.pred$class)
##           predicted
## correct    low med_low med_high high
##   low       14       9        1    0
##   med_low    2      21        5    0
##   med_high   0       9       20    0
##   high       0       0        1   20

We saved the crime categories from the test set and then removed the categorical crime variable from the test dataset. Then we predicted the classes with the LDA model on the test data. We also cross tabulated the results from prediction with the crime categories from the test set. We notice that the predicted values of the class variable that are predicted on the test data using the results from the LDA on the train data differ partly from the observed values of the class variable in the test data. Anyhow the predictions are same as the observed values more frequently that the predictions differ from observed values. Prediction is hence quite good.

# load the data
data("Boston")

# center and standardize variables
boston_scaled <- scale(Boston)

# change the object to data frame
boston_scaled=as.data.frame(boston_scaled)

# euclidean distance matrix
dist_eu <- dist(boston_scaled)

# look at the summary of the distances
summary(dist_eu)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1343  3.4620  4.8240  4.9110  6.1860 14.4000
# manhattan distance matrix
dist_man <- dist(boston_scaled, method="manhattan")

# look at the summary of the distances

summary(dist_man)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2662  8.4830 12.6100 13.5500 17.7600 48.8600
# k-means clustering
km <-kmeans(dist_eu, centers = 15)

# plot the Boston dataset with clusters
pairs(boston_scaled, col = km$cluster)

set.seed(123)


# determine the number of clusters
k_max <- 15

# calculate the total within sum of squares
twcss <- sapply(1:k_max, function(k){kmeans(dist_eu, k)$tot.withinss})

# visualize the results
plot(1:k_max, twcss, type='b')

# k-means clustering
km <-kmeans(dist_eu, centers = 2)

# plot the Boston dataset with clusters
pairs(boston_scaled, col = km$cluster)

We reloaded the Boston dataset and standardized the dataset. We calculated the distances between the observations using the euclidian and manhattan distance matrixes. We ran k-means algorithm on the dataset and visualized the clusters with pairs-function where clusters are separated with colors. Then we investigated what is the optimal number of clusters is with total within sum of squares. We drew a picture where the number of clusters in on the x-axis and quantity of total within sum of squares on the y-axis. From the picture we saw that total sum of squares gets its highest value at one. We thought one is too small number of clusters and chose two for the number of clusters. We ran the k-means algorithm again and chose two as the number of centers. We visualized the clusters again with pairs-function.

Super-bonus-exercise

model_predictors <- dplyr::select(train, -crime)

# check the dimensions
dim(model_predictors)
## [1] 404  13
dim(lda.fit$scaling)
## [1] 13  3
# matrix multiplication
matrix_product <- as.matrix(model_predictors) %*% lda.fit$scaling
matrix_product <- as.data.frame(matrix_product)

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
plot_ly(x = matrix_product$LD1, y = matrix_product$LD2, z = matrix_product$LD3, type= 'scatter3d', mode='markers')
plot_ly(x = matrix_product$LD1, y = matrix_product$LD2, z = matrix_product$LD3, type= 'scatter3d', mode='markers', color=train$crime)

We created a matrix product, which is a projection of the data points. We created a 3D plot of the columns of the matrix product. We added argument color as a argument in the plot_ly() function and set the color to be the crime classes of the train set. Hence we got two plots. In the second plot, crime classes low, medium low, medium high and high are separated with colors. Otherwise two plots are similar.